Data Visualization

Questions

  • How do I visualize results using open source software?

Objectives

  • Learn about Holoviews for data visualization

import numpy as np
import xarray as xr
import holoviews as hv
from holoviews import opts
import geoviews as gv
import datashader as ds
import cartopy.crs as ccrs
import xarray as xa
from holoviews.operation.datashader import regrid, shade
from bokeh.tile_providers import OSM

%matplotlib inline
hv.extension('bokeh', width=80)
---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
<ipython-input-1-48edf7250545> in <module>
----> 1 import numpy as np
      2 import xarray as xr
      3 import holoviews as hv
      4 from holoviews import opts
      5 import geoviews as gv

ModuleNotFoundError: No module named 'numpy'

Check for single band

url = 'http://landsat-pds.s3.amazonaws.com/c1/L8/227/065/LC08_L1TP_227065_20200608_20200626_01_T1/'
redband = url+'LC08_L1TP_227065_20200608_20200626_01_T1_B{}.TIF'.format(4)

redband
'http://landsat-pds.s3.amazonaws.com/c1/L8/227/065/LC08_L1TP_227065_20200608_20200626_01_T1/LC08_L1TP_227065_20200608_20200626_01_T1_B4.TIF'

Read in all the bands

file_path = 'http://landsat-pds.s3.amazonaws.com/c1/L8/227/065/LC08_L1TP_227065_20200608_20200626_01_T1/LC08_L1TP_227065_20200608_20200626_01_T1_B%s.TIF'
bands = list(range(1, 12))
bands = [xr.open_rasterio(file_path%band).load()[0] for band in bands]
bands
[<xarray.DataArray (y: 7761, x: 7621)>
 array([[0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        ...,
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0]], dtype=uint16)
 Coordinates:
     band     int64 1
   * y        (y) float64 -6.837e+05 -6.837e+05 ... -9.165e+05 -9.165e+05
   * x        (x) float64 5.733e+05 5.733e+05 5.734e+05 ... 8.019e+05 8.019e+05
 Attributes:
     transform:      (30.0, 0.0, 573285.0, 0.0, -30.0, -683685.0)
     crs:            +init=epsg:32621
     res:            (30.0, 30.0)
     is_tiled:       1
     nodatavals:     (nan,)
     scales:         (1.0,)
     offsets:        (0.0,)
     AREA_OR_POINT:  Point,
 <xarray.DataArray (y: 7761, x: 7621)>
 array([[0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        ...,
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0]], dtype=uint16)
 Coordinates:
     band     int64 1
   * y        (y) float64 -6.837e+05 -6.837e+05 ... -9.165e+05 -9.165e+05
   * x        (x) float64 5.733e+05 5.733e+05 5.734e+05 ... 8.019e+05 8.019e+05
 Attributes:
     transform:      (30.0, 0.0, 573285.0, 0.0, -30.0, -683685.0)
     crs:            +init=epsg:32621
     res:            (30.0, 30.0)
     is_tiled:       1
     nodatavals:     (nan,)
     scales:         (1.0,)
     offsets:        (0.0,)
     AREA_OR_POINT:  Point,
 <xarray.DataArray (y: 7761, x: 7621)>
 array([[0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        ...,
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0]], dtype=uint16)
 Coordinates:
     band     int64 1
   * y        (y) float64 -6.837e+05 -6.837e+05 ... -9.165e+05 -9.165e+05
   * x        (x) float64 5.733e+05 5.733e+05 5.734e+05 ... 8.019e+05 8.019e+05
 Attributes:
     transform:      (30.0, 0.0, 573285.0, 0.0, -30.0, -683685.0)
     crs:            +init=epsg:32621
     res:            (30.0, 30.0)
     is_tiled:       1
     nodatavals:     (nan,)
     scales:         (1.0,)
     offsets:        (0.0,)
     AREA_OR_POINT:  Point,
 <xarray.DataArray (y: 7761, x: 7621)>
 array([[0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        ...,
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0]], dtype=uint16)
 Coordinates:
     band     int64 1
   * y        (y) float64 -6.837e+05 -6.837e+05 ... -9.165e+05 -9.165e+05
   * x        (x) float64 5.733e+05 5.733e+05 5.734e+05 ... 8.019e+05 8.019e+05
 Attributes:
     transform:      (30.0, 0.0, 573285.0, 0.0, -30.0, -683685.0)
     crs:            +init=epsg:32621
     res:            (30.0, 30.0)
     is_tiled:       1
     nodatavals:     (nan,)
     scales:         (1.0,)
     offsets:        (0.0,)
     AREA_OR_POINT:  Point,
 <xarray.DataArray (y: 7761, x: 7621)>
 array([[0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        ...,
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0]], dtype=uint16)
 Coordinates:
     band     int64 1
   * y        (y) float64 -6.837e+05 -6.837e+05 ... -9.165e+05 -9.165e+05
   * x        (x) float64 5.733e+05 5.733e+05 5.734e+05 ... 8.019e+05 8.019e+05
 Attributes:
     transform:      (30.0, 0.0, 573285.0, 0.0, -30.0, -683685.0)
     crs:            +init=epsg:32621
     res:            (30.0, 30.0)
     is_tiled:       1
     nodatavals:     (nan,)
     scales:         (1.0,)
     offsets:        (0.0,)
     AREA_OR_POINT:  Point,
 <xarray.DataArray (y: 7761, x: 7621)>
 array([[0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        ...,
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0]], dtype=uint16)
 Coordinates:
     band     int64 1
   * y        (y) float64 -6.837e+05 -6.837e+05 ... -9.165e+05 -9.165e+05
   * x        (x) float64 5.733e+05 5.733e+05 5.734e+05 ... 8.019e+05 8.019e+05
 Attributes:
     transform:      (30.0, 0.0, 573285.0, 0.0, -30.0, -683685.0)
     crs:            +init=epsg:32621
     res:            (30.0, 30.0)
     is_tiled:       1
     nodatavals:     (nan,)
     scales:         (1.0,)
     offsets:        (0.0,)
     AREA_OR_POINT:  Point,
 <xarray.DataArray (y: 7761, x: 7621)>
 array([[0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        ...,
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0]], dtype=uint16)
 Coordinates:
     band     int64 1
   * y        (y) float64 -6.837e+05 -6.837e+05 ... -9.165e+05 -9.165e+05
   * x        (x) float64 5.733e+05 5.733e+05 5.734e+05 ... 8.019e+05 8.019e+05
 Attributes:
     transform:      (30.0, 0.0, 573285.0, 0.0, -30.0, -683685.0)
     crs:            +init=epsg:32621
     res:            (30.0, 30.0)
     is_tiled:       1
     nodatavals:     (nan,)
     scales:         (1.0,)
     offsets:        (0.0,)
     AREA_OR_POINT:  Point,
 <xarray.DataArray (y: 15521, x: 15241)>
 array([[0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        ...,
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0]], dtype=uint16)
 Coordinates:
     band     int64 1
   * y        (y) float64 -6.837e+05 -6.837e+05 ... -9.165e+05 -9.165e+05
   * x        (x) float64 5.733e+05 5.733e+05 5.733e+05 ... 8.019e+05 8.019e+05
 Attributes:
     transform:      (15.0, 0.0, 573292.5, 0.0, -15.0, -683692.5)
     crs:            +init=epsg:32621
     res:            (15.0, 15.0)
     is_tiled:       1
     nodatavals:     (nan,)
     scales:         (1.0,)
     offsets:        (0.0,)
     AREA_OR_POINT:  Point,
 <xarray.DataArray (y: 7761, x: 7621)>
 array([[0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        ...,
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0]], dtype=uint16)
 Coordinates:
     band     int64 1
   * y        (y) float64 -6.837e+05 -6.837e+05 ... -9.165e+05 -9.165e+05
   * x        (x) float64 5.733e+05 5.733e+05 5.734e+05 ... 8.019e+05 8.019e+05
 Attributes:
     transform:      (30.0, 0.0, 573285.0, 0.0, -30.0, -683685.0)
     crs:            +init=epsg:32621
     res:            (30.0, 30.0)
     is_tiled:       1
     nodatavals:     (nan,)
     scales:         (1.0,)
     offsets:        (0.0,)
     AREA_OR_POINT:  Point,
 <xarray.DataArray (y: 7761, x: 7621)>
 array([[0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        ...,
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0]], dtype=uint16)
 Coordinates:
     band     int64 1
   * y        (y) float64 -6.837e+05 -6.837e+05 ... -9.165e+05 -9.165e+05
   * x        (x) float64 5.733e+05 5.733e+05 5.734e+05 ... 8.019e+05 8.019e+05
 Attributes:
     transform:      (30.0, 0.0, 573285.0, 0.0, -30.0, -683685.0)
     crs:            +init=epsg:32621
     res:            (30.0, 30.0)
     is_tiled:       1
     nodatavals:     (nan,)
     scales:         (1.0,)
     offsets:        (0.0,)
     AREA_OR_POINT:  Point,
 <xarray.DataArray (y: 7761, x: 7621)>
 array([[0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        ...,
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0]], dtype=uint16)
 Coordinates:
     band     int64 1
   * y        (y) float64 -6.837e+05 -6.837e+05 ... -9.165e+05 -9.165e+05
   * x        (x) float64 5.733e+05 5.733e+05 5.734e+05 ... 8.019e+05 8.019e+05
 Attributes:
     transform:      (30.0, 0.0, 573285.0, 0.0, -30.0, -683685.0)
     crs:            +init=epsg:32621
     res:            (30.0, 30.0)
     is_tiled:       1
     nodatavals:     (nan,)
     scales:         (1.0,)
     offsets:        (0.0,)
     AREA_OR_POINT:  Point]
opts.defaults(opts.RGB(width=600, height=600))

nodata= 1

def one_band(b):
    xs, ys = b['x'], b['y']
    b = ds.utils.orient_array(b)
    a = (np.where(np.logical_or(np.isnan(b),b<=nodata),0,255)).astype(np.uint8)    
    col, rows = b.shape
    return hv.RGB((xs, ys[::-1], b, b, b, a), vdims=list('RGBA'))

# not working through GeoViews
#tiles = gv.WMTS(OSM)
tiles = hv.element.tiles.EsriUSATopo().opts(width=600, height=550)
tiles*shade(regrid(one_band(bands[1])), cmap=['black', 'white']).redim(x='Longitude', y='Latitude')

Add background tiles

from datashader.utils import ngjit
nodata= 1

@ngjit
def normalize_data(agg):
    out = np.zeros_like(agg)
    min_val = 0
    max_val = 2**16 - 1
    range_val = max_val - min_val
    col, rows = agg.shape
    c = 40
    th = .125
    for x in range(col):
        for y in range(rows):
            val = agg[x, y]
            norm = (val - min_val) / range_val
            norm = 1 / (1 + np.exp(c * (th - norm))) # bonus
            out[x, y] = norm * 255.0
    return out

def combine_bands(r, g, b):
    xs, ys = r['x'], r['y']
    r, g, b = [ds.utils.orient_array(img) for img in (r, g, b)]
    a = (np.where(np.logical_or(np.isnan(r),r<=nodata),0,255)).astype(np.uint8)    
    r = (normalize_data(r)).astype(np.uint8)
    g = (normalize_data(g)).astype(np.uint8)
    b = (normalize_data(b)).astype(np.uint8)
    return hv.RGB((xs, ys[::-1], r, g, b, a), vdims=list('RGBA'))
# Lets visualize area of in true color
true_color = combine_bands(bands[3], bands[2], bands[1]).relabel("True Color (R=Red, G=Green, B=Blue)")
regrid(true_color)
import pandas as pd
band_info = pd.DataFrame([
        (1,  "Aerosol", " 0.43 - 0.45",    0.440,  "30",         "Coastal aerosol"),
        (2,  "Blue",    " 0.45 - 0.51",    0.480,  "30",         "Blue"),
        (3,  "Green",   " 0.53 - 0.59",    0.560,  "30",         "Green"),
        (4,  "Red",     " 0.64 - 0.67",    0.655,  "30",         "Red"),
        (5,  "NIR",     " 0.85 - 0.88",    0.865,  "30",         "Near Infrared (NIR)"),
        (6,  "SWIR1",   " 1.57 - 1.65",    1.610,  "30",         "Shortwave Infrared (SWIR) 1"),
        (7,  "SWIR2",   " 2.11 - 2.29",    2.200,  "30",         "Shortwave Infrared (SWIR) 2"),
        (8,  "Panc",    " 0.50 - 0.68",    0.590,  "15",         "Panchromatic"),
        (9,  "Cirrus",  " 1.36 - 1.38",    1.370,  "30",         "Cirrus"),
        (10, "TIRS1",   "10.60 - 11.19",   10.895, "100 * (30)", "Thermal Infrared (TIRS) 1"),
        (11, "TIRS2",   "11.50 - 12.51",   12.005, "100 * (30)", "Thermal Infrared (TIRS) 2")],
    columns=['Band', 'Name', 'Wavelength Range (µm)', 'Nominal Wavelength (µm)', 'Resolution (m)', 'Description']).set_index(["Band"])
band_info
Name Wavelength Range (µm) Nominal Wavelength (µm) Resolution (m) Description
Band
1 Aerosol 0.43 - 0.45 0.440 30 Coastal aerosol
2 Blue 0.45 - 0.51 0.480 30 Blue
3 Green 0.53 - 0.59 0.560 30 Green
4 Red 0.64 - 0.67 0.655 30 Red
5 NIR 0.85 - 0.88 0.865 30 Near Infrared (NIR)
6 SWIR1 1.57 - 1.65 1.610 30 Shortwave Infrared (SWIR) 1
7 SWIR2 2.11 - 2.29 2.200 30 Shortwave Infrared (SWIR) 2
8 Panc 0.50 - 0.68 0.590 15 Panchromatic
9 Cirrus 1.36 - 1.38 1.370 30 Cirrus
10 TIRS1 10.60 - 11.19 10.895 100 * (30) Thermal Infrared (TIRS) 1
11 TIRS2 11.50 - 12.51 12.005 100 * (30) Thermal Infrared (TIRS) 2
combos = pd.DataFrame([
    (4,3,2,"True color",""),
    (7,6,4,"Urban","False color"),
    (5,4,3,"Vegetation","Color Infrared"),
    (6,5,2,"Agriculture",""),
    (7,6,5,"Penetration","Atmospheric Penetration"),
    (5,6,2,"Healthy Vegetation",""),
    (5,6,4,"Land vs. Water",""),
    (7,5,3,"Atmosphere Removal","Natural With Atmospheric Removal"),
    (7,5,4,"Shortwave Infrared",""),
    (6,5,4,"Vegetation Analysis","")],
    columns=['R', 'G', 'B', 'Name', 'Description']).set_index(["Name"])
combos

def combo(name):
    c=combos.loc[name]
    return regrid(combine_bands(bands[c.R-1],bands[c.G-1],bands[c.B-1])).relabel(name)

layout = combo("Urban") + combo("Vegetation") + combo("Agriculture") + combo("Land vs. Water")

layout.opts(
    opts.RGB(width=350, height=350, xaxis=None, yaxis=None, framewise=True)).cols(2)

Explore the full hyperspectral information at any location (here Amazon forest) in the true-color image. Click the point on the left- updates the spectral curve whenever you hover over an area of the image:

band_map = hv.HoloMap({i: hv.Image(band) for i, band in enumerate(bands)})

def spectrum(x, y):
    try: 
        spectrum_vals = band_map.sample(x=x, y=y)['z'][:-1]
        point = gv.Points([(x, y)], crs=ccrs.GOOGLE_MERCATOR)
        point = gv.operation.project_points(point, projection=ccrs.PlateCarree())
        label = label = 'Lon: %.3f, Lat: %.3f' % tuple(point.array()[0])
    except:
        spectrum_vals = np.zeros(11)
        label = 'Lon: -, Lat: -'
    
    return hv.Curve((band_info['Nominal Wavelength (µm)'].values, spectrum_vals), label=label,
                    kdims=['Wavelength (µm)'], vdims=['Luminance']).sort()

# x and y give the location in Web Mercator coordinates
spectrum(x=700000, y=-800000).opts(width=800, height=300, logx=True)
WARNING:param.HoloMap05936: The HoloMap.sample method is deprecated, for equivalent functionality use HoloMap.apply.sample().collapse().
tap = hv.streams.PointerXY(source=true_color)
spectrum_curve = hv.DynamicMap(spectrum, streams=[tap]).redim.range(Luminance=(0, 30000))

layout = regrid(true_color) + spectrum_curve
layout.opts(
    opts.Curve(width=450, height=450, logx=True),
    opts.RGB(width=450, height=450))
WARNING:param.HoloMap05936: The HoloMap.sample method is deprecated, for equivalent functionality use HoloMap.apply.sample().collapse().
# Can change the background tiles now
# see here https://holoviews.org/reference/elements/bokeh/Tiles.html
tiles = hv.element.tiles.StamenWatercolor().opts(width=600, height=550)
layout = tiles * regrid(true_color) + spectrum_curve
layout.opts(
    opts.Curve(width=450, height=450, logx=True),
    opts.RGB(width=450, height=450))

Point it to particular location

Check ESRI imagery

hv.element.tiles.EsriImagery().opts(width=600, height=550)
#lets visualize the ndvi
ndvi = (bands[5]-bands[4]) / (bands[5] + bands[4])
def ndvi(r, nir):  
    r = (normalize_data(r))
    nir = (normalize_data(nir))
    return((r-nir)/(r+nir))